/**
 * @file AM_RAYFAC.cpp
 * @brief Computes curvature radii of contact facets
 * 
 * This module calculates the curvature radii for contact facets in adaptive mesh refinement.
 * 
 */

#include <vector>
#include <algorithm>
#include <cmath>
#include <memory>

// Type definitions
using Real = float;      // Single precision floating point
using Integer = int;     // Integer type
using RealArray = std::vector<Real>;
using IntArray = std::vector<Integer>;

// Constants
constexpr Real RRMAX = 1.0e31f;  // Maximum radius value
constexpr Integer ITOW = 1;       // Task increment factor

// Global variables (would normally be in a namespace or class)
Integer NFA;            // Number of facets
Integer IFACAD;         // Adaptive facets flag
Integer IDMP;           // DMP flag
Integer MPROC;          // Number of processors
Integer NPROC;          // Number of processes
Integer MXTASK;         // Maximum number of tasks
Integer ICONTA;         // Contact flag
Integer ICTTHK;         // Thickness flag
Integer KTYPF4;         // Facet type 4 flag

// Memory management class
class MemoryManager {
public:
    static void allocate(const std::string& name, size_t size, int type) {
        // Implementation would handle memory allocation
    }
    
    static void free(const std::string& name) {
        // Implementation would handle memory deallocation
    }
    
    static void resize(const std::string& name, size_t newSize, int type) {
        // Implementation would handle memory resizing
    }
};

// Main function to compute curvature radii
void AM_RAYFAC(RealArray& MA) {
    // =====================================================================
    // Build edges array of contact facets
    // =====================================================================
    Integer NTEMP = 4 * NFA;
    Integer NEDGES = 0;
    Integer NFATOT = 0;
    
    // Allocate temporary arrays for edge processing
    MemoryManager::allocate("KT_INO1", NTEMP, 2);
    MemoryManager::allocate("KT_INO2", NTEMP, 2);
    MemoryManager::allocate("KT_ITRART", NTEMP, 2);
    MemoryManager::allocate("KT_IADRES", NTEMP, 2);
    MemoryManager::allocate("KT_IFLART", NTEMP, 2);
    MemoryManager::allocate("KT_NELART", NTEMP, 2);
    MemoryManager::allocate("KT_INOAR1", NTEMP, 2);
    MemoryManager::allocate("KT_INOAR2", NTEMP, 2);
    MemoryManager::allocate("KT_IARPER", NTEMP, 2);
    MemoryManager::allocate("KT_IART12", NFA, 2);
    MemoryManager::allocate("KT_IART23", NFA, 2);
    MemoryManager::allocate("KT_IART34", NFA, 2);
    MemoryManager::allocate("KT_IART41", NFA, 2);
    
    // Build the edge arrays
    AM_MAKART(
        MA[LB("KT_INO1")], MA[LB("KT_INO2")],
        MA[LB("KT_ITRART")], MA[LB("KT_IADRES")],
        MA[LB("KT_IFLART")], MA[LB("KT_NELART")],
        MA[LB("KT_INOAR1")], MA[LB("KT_INOAR2")],
        MA[LB("KT_IARPER")], MA[LB("KT_IART12")],
        MA[LB("KT_IART23")], MA[LB("KT_IART34")],
        MA[LB("KT_IART41")], MA[LB("KI_ICFA1")],
        MA[LB("KI_ICFA2")], MA[LB("KI_ICFA3")],
        MA[LB("KI_ICFA4")], MA[LB("KI_NFINUM")],
        MA[LB("KI_NFITYP")], MA[LB("KI_NFIPT")],
        MA[LB("KI_NFINPR")], MA[LB("KI_NFIPER")],
        NEDGES, NFATOT
    );
    
    // Free temporary arrays
    MemoryManager::free("KT_INO1");
    MemoryManager::free("KT_INO2");
    MemoryManager::free("KT_ITRART");
    MemoryManager::free("KT_IADRES");
    MemoryManager::free("KT_IFLART");
    
    // =====================================================================
    // Compute curvature radius on contact facets
    // =====================================================================
    // Free old curvature arrays
    MemoryManager::free("KR_RCOURB");
    MemoryManager::free("KR_DISTAR");
    
    // Resize edge arrays
    MemoryManager::resize("KT_INOAR1", NEDGES, 2);
    MemoryManager::resize("KT_INOAR2", NEDGES, 2);
    MemoryManager::resize("KT_IARPER", NEDGES, 2);
    
    // MPP: build contact edge mailing using ID nodes
    if (IDMP == 1) {
        MemoryManager::allocate("KT_IDNOAR1", NEDGES, 2);
        MemoryManager::allocate("KT_IDNOAR2", NEDGES, 2);
        MemoryManager::allocate("KT_ITRART", NEDGES, 2);
        
        MPP_MMAILCTEDG(
            MA[LB("KT_ITRART")], MA[LB("KT_INOAR1")], MA[LB("KT_INOAR2")],
            MA[LB("KI_IDN1CTEDA")], MA[LB("KI_IDN2CTEDA")], MA[LB("KI_IDN1CTEDB")],
            MA[LB("KI_IDN2CTEDB")], MA[LB("KT_IDNOAR1")], MA[LB("KT_IDNOAR2")],
            MA[LB("KI_MAILCTEDA")], MA[LB("KI_MAILCTEDB")], MA[LB("KI_IDNOD")],
            NEDGES
        );
        
        MemoryManager::free("KT_IDNOAR1");
        MemoryManager::free("KT_IDNOAR2");
        MemoryManager::free("KT_ITRART");
    }
    
    // Allocate arrays for facet processing
    MemoryManager::allocate("KT_IPTART", NEDGES+1, 2);
    MemoryManager::allocate("KT_IFACAR", NFATOT, 2);
    MemoryManager::allocate("KT_ITYFAC", NFATOT, 2);
    MemoryManager::allocate("KT_XGFA", NFA, 3);
    MemoryManager::allocate("KT_YGFA", NFA, 3);
    MemoryManager::allocate("KT_ZGFA", NFA, 3);
    MemoryManager::allocate("KT_VNX", NFA, 3);
    MemoryManager::allocate("KT_VNY", NFA, 3);
    MemoryManager::allocate("KT_VNZ", NFA, 3);
    MemoryManager::allocate("KT_ITYCFA", NFA, 2);
    MemoryManager::allocate("KT_IDELFA", NFA, 2);
    
    // Copy facet type information
    std::copy_n(
        MA.begin() + LB("KI_ITYCFA"), 
        NFA, 
        MA.begin() + LB("KT_ITYCFA")
    );
    
    // Process edges and compute normals
    AM_RGARET1(
        MA[LB("KT_NELART")], MA[LB("KT_IARPER")],
        MA[LB("KT_IART12")], MA[LB("KT_IART23")],
        MA[LB("KT_IART34")], MA[LB("KT_IART41")],
        MA[LB("KT_IPTART")], MA[LB("KT_IFACAR")],
        MA[LB("KT_ITYFAC")], MA[LB("KT_VNX")],
        MA[LB("KT_VNY")], MA[LB("KT_VNZ")],
        MA[LB("KI_ICFA1")], MA[LB("KI_ICFA2")],
        MA[LB("KI_ICFA3")], MA[LB("KI_ICFA4")],
        MA[LB("KR_X")], MA[LB("KT_XGFA")],
        MA[LB("KT_YGFA")], MA[LB("KT_ZGFA")],
        MA[LB("KT_ITYCFA")], MA[LB("KT_IDELFA")],
        MA[LB("KI_IDSHE")], MA[LB("KI_IDBRI")],
        MA[LB("KI_IELCFA")], NEDGES
    );
    
    // MPP: Communication of NELART from other procs
    if (IDMP == 1) {
        MemoryManager::allocate("KT_NELARTADD", NEDGES, 2);
        MemoryManager::allocate("KT_IPTARTADD", LGMAILCTEDA+1, 2);
        MemoryManager::allocate("KT_NFAPROCR", MPROC, 2);
        MemoryManager::allocate("KT_IPTFAPRO", MPROC+1, 2);
        MemoryManager::allocate("KT_IFALOCFA", NFA, 2);
        
        Integer LIBUFMPI = std::max(LGMAILCTEDA, LGMAILCTEDB) + 1;
        MemoryManager::allocate("KT_IBUFMPI", LIBUFMPI, 2);
        
        // Count number of information to transfer
        Integer NFAADD = 0;
        Integer NFAADDS = 0;
        Integer NFAADDTOT = 0;
        Integer NFAADDTOTS = 0;
        
        COMM_CTEDG_LEN(
            MA[LB("KT_NELART")], MA[LB("KT_NELARTADD")], MA[LB("KT_IPTART")],
            MA[LB("KI_LCTEDBPA")], MA[LB("KI_IPTCTEDRA")], MA[LB("KI_MAILCTEDA")],
            MA[LB("KI_LCTEDBPB")], MA[LB("KI_IPTCTEDRB")], MA[LB("KI_MAILCTEDB")],
            MA[LB("KT_IBUFMPI")], MA[LB("KT_IFACAR")], MA[LB("KT_NFAPROCR")],
            MA[LB("KT_IFALOCFA")], MA[LB("KT_IPTFAPRO")], MA[LB("KT_IPTARTADD")],
            NFAADD, NFAADDTOT, NFAADDTOTS, NFAADDS
        );
        
        MemoryManager::free("KT_IBUFMPI");
        
        // Allocate arrays for additional facets
        MemoryManager::allocate("KT_VNXADD", NFAADD, 3);
        MemoryManager::allocate("KT_VNYADD", NFAADD, 3);
        MemoryManager::allocate("KT_VNZADD", NFAADD, 3);
        MemoryManager::allocate("KT_XGFAADD", NFAADD, 3);
        MemoryManager::allocate("KT_YGFAADD", NFAADD, 3);
        MemoryManager::allocate("KT_ZGFAADD", NFAADD, 3);
        MemoryManager::allocate("KT_ITYFACADD", NFAADD, 2);
        MemoryManager::allocate("KT_ITYCFAADD", NFAADD, 2);
        MemoryManager::allocate("KT_IDELFAADD", NFAADD, 2);
        MemoryManager::allocate("KT_IFAFALOC", NFA, 2);
        
        MemoryManager::allocate("KT_IFACARADD", NFAADDTOT, 2);
        LIBUFMPI = std::max(
            NFAADDTOT + (6*ITOW+3)*NFAADD,
            NFAADDTOTS + (6*ITOW+3)*NFAADDS
        );
        MemoryManager::allocate("KT_IBUFMPI", LIBUFMPI, 2);
        
        // Send information between processors
        COMM_CTEDG_VAR(
            MA[LB("KT_NELART")], MA[LB("KT_NELARTADD")], MA[LB("KT_IFACAR")],
            MA[LB("KT_IPTART")], MA[LB("KT_IPTARTADD")], MA[LB("KT_IBUFMPI")],
            MA[LB("KT_VNX")], MA[LB("KT_VNY")], MA[LB("KT_VNZ")],
            MA[LB("KT_VNXADD")], MA[LB("KT_VNYADD")], MA[LB("KT_VNZADD")],
            MA[LB("KT_XGFA")], MA[LB("KT_YGFA")], MA[LB("KT_ZGFA")],
            MA[LB("KT_XGFAADD")], MA[LB("KT_YGFAADD")], MA[LB("KT_ZGFAADD")],
            MA[LB("KT_ITYFAC")], MA[LB("KI_ITYCFA")], MA[LB("KI_IELCFA")],
            MA[LB("KI_IDSHE")], MA[LB("KI_IDBRI")],
            MA[LB("KT_ITYFACADD")], MA[LB("KT_ITYCFAADD")], MA[LB("KT_IDELFAADD")],
            MA[LB("KI_LCTEDBPA")], MA[LB("KI_IPTCTEDRA")], MA[LB("KI_MAILCTEDA")],
            MA[LB("KI_LCTEDBPB")], MA[LB("KI_IPTCTEDRB")], MA[LB("KI_MAILCTEDB")],
            MA[LB("KT_IFALOCFA")], MA[LB("KT_IFAFALOC")], MA[LB("KT_IPTFAPRO")],
            MA[LB("KT_NFAPROCR")], MA[LB("KT_IFACARADD")], LIBUFMPI
        );
        
        MemoryManager::free("KT_IBUFMPI");
        
        // Resize and update pointers after adding new information
        Integer LNEW = NFATOT + NFAADDTOT;
        MemoryManager::resize("KT_IFACAR", LNEW, 2);
        LNEW = NFA + NFAADD;
        MemoryManager::resize("KT_ITYFAC", LNEW, 2);
        MemoryManager::resize("KT_XGFA", LNEW, 3);
        MemoryManager::resize("KT_YGFA", LNEW, 3);
        MemoryManager::resize("KT_ZGFA", LNEW, 3);
        MemoryManager::resize("KT_VNX", LNEW, 3);
        MemoryManager::resize("KT_VNY", LNEW, 3);
        MemoryManager::resize("KT_VNZ", LNEW, 3);
        MemoryManager::resize("KT_ITYCFA", LNEW, 2);
        MemoryManager::resize("KT_IDELFA", LNEW, 2);
        
        // Save old pointers
        MemoryManager::allocate("KT_IPTARTOLD", NEDGES+1, 2);
        MemoryManager::allocate("KT_IFACAROLD", NFATOT, 2);
        std::copy_n(
            MA.begin() + LB("KT_IPTART"), 
            NEDGES+1, 
            MA.begin() + LB("KT_IPTARTOLD")
        );
        std::copy_n(
            MA.begin() + LB("KT_IFACAR"), 
            NFATOT, 
            MA.begin() + LB("KT_IFACAROLD")
        );
        
        // Update with new information
        AM_RGARET2(
            MA[LB("KT_NELART")], MA[LB("KT_IPTART")], MA[LB("KT_IFACAR")],
            MA[LB("KI_MAILCTEDA")], MA[LB("KT_NELARTADD")], MA[LB("KT_IPTARTADD")],
            MA[LB("KT_IFACARADD")], MA[LB("KT_ITYFAC")], MA[LB("KT_VNX")],
            MA[LB("KT_VNY")], MA[LB("KT_VNZ")], MA[LB("KT_ITYFACADD")],
            MA[LB("KT_VNXADD")], MA[LB("KT_VNYADD")], MA[LB("KT_VNZADD")],
            MA[LB("KT_XGFA")], MA[LB("KT_YGFA")], MA[LB("KT_ZGFA")],
            MA[LB("KT_XGFAADD")], MA[LB("KT_YGFAADD")], MA[LB("KT_ZGFAADD")],
            MA[LB("KT_IDELFA")], MA[LB("KT_IDELFAADD")], MA[LB("KT_IPTARTOLD")],
            MA[LB("KT_IFACAROLD")], NEDGES
        );
        
        // Free temporary arrays
        MemoryManager::free("KT_NFAPROCR");
        MemoryManager::free("KT_IFALOCFA");
        MemoryManager::free("KT_IFAFALOC");
        MemoryManager::free("KT_IPTFAPRO");
        MemoryManager::free("KT_IPTARTOLD");
        MemoryManager::free("KT_IFACAROLD");
    }
    
    // =====================================================================
    // Compute curvature radii
    // =====================================================================
    MemoryManager::allocate("KR_RCOURB", NEDGES, 3);
    MemoryManager::allocate("KR_DISTAR", NEDGES, 3);
    
    // Initialize with maximum values
    std::fill_n(MA.begin() + LB("KR_RCOURB"), NEDGES, RRMAX);
    std::fill_n(MA.begin() + LB("KR_DISTAR"), NEDGES, RRMAX);
    
    // Parallel computation of curvature radii
    Integer LGRAIN = LGR; // Grain size for parallel processing
    Integer NTASK = 0;
    Integer NIL = 0;
    Integer ISL = 0;
    
    // Distribute work among tasks (simplified version)
    DISPAR(NIL, ISL, NEDGES, NPROC, LGRAIN, LGRAIN, MXTASK, NTASK);
    
    // Parallel loop to compute curvature radii
    #pragma omp parallel for
    for (Integer I = 0; I < NTASK; ++I) {
        Integer JSL = ISL * ITOW;
        AM_CALRAY(
            MA[LB("KT_INOAR1") + ISL], MA[LB("KT_INOAR2") + ISL],
            MA[LB("KT_IPTART") + ISL], MA[LB("KT_IFACAR")],
            MA[LB("KT_ITYFAC")], MA[LB("KT_VNX")],
            MA[LB("KT_VNY")], MA[LB("KT_VNZ")],
            MA[LB("KR_RCOURB") + JSL], MA[LB("KR_DISTAR") + JSL],
            MA[LB("KI_ICFA1")], MA[LB("KI_ICFA2")],
            MA[LB("KI_ICFA3")], MA[LB("KI_ICFA4")],
            MA[LB("KT_ITYCFA")], MA[LB("KT_IDELFA")],
            MA[LB("KR_X")], MA[LB("KR_THKCMP")],
            MA[LB("KT_XGFA")], MA[LB("KT_YGFA")],
            MA[LB("KT_ZGFA")], NIL,
            KTYPF4, ICONTA, ICTTHK
        );
    }
    
    // Free temporary arrays
    MemoryManager::free("KT_IPTART");
    MemoryManager::free("KT_IFACAR");
    MemoryManager::free("KT_ITYFAC");
    MemoryManager::free("KT_VNX");
    MemoryManager::free("KT_VNY");
    MemoryManager::free("KT_VNZ");
    MemoryManager::free("KT_XGFA");
    MemoryManager::free("KT_YGFA");
    MemoryManager::free("KT_ZGFA");
    MemoryManager::free("KT_ITYCFA");
    MemoryManager::free("KT_IDELFA");
    
    if (IDMP == 1) {
        MemoryManager::free("KT_NELARTADD");
        MemoryManager::free("KT_IPTARTADD");
        MemoryManager::free("KT_VNXADD");
        MemoryManager::free("KT_VNYADD");
        MemoryManager::free("KT_VNZADD");
        MemoryManager::free("KT_XGFAADD");
        MemoryManager::free("KT_YGFAADD");
        MemoryManager::free("KT_ZGFAADD");
        MemoryManager::free("KT_ITYFACADD");
        MemoryManager::free("KT_ITYCFAADD");
        MemoryManager::free("KT_IDELFAADD");
        MemoryManager::free("KT_IFACARADD");
    }
    
    // MPP: Send back curvature values
    if (IDMP == 1) {
        Integer LRBUFMPI = 2 * std::max(LGMAILCTEDA, LGMAILCTEDB);
        MemoryManager::allocate("KT_RBUFMPI", LRBUFMPI, 3);
        Integer ITAG = 2442;
        
        COMM_RP2P(
            ITAG, MA[LB("KT_RBUFMPI")], LRBUFMPI,
            NCTEDBPA, MA[LB("KI_LCTEDBPA")], MA[LB("KI_IPTCTEDRA")],
            MA[LB("KI_MAILCTEDA")],
            NCTEDBPB, MA[LB("KI_LCTEDBPB")], MA[LB("KI_IPTCTEDRB")],
            MA[LB("KI_MAILCTEDB")],
            MA[LB("KR_RCOURB")], 1, 1, "SENDRECV",
            MA[LB("KR_DISTAR")], 1, 1, "SENDRECV"
        );
        
        MemoryManager::free("KT_RBUFMPI");
    }
    
    // Modification on father's edges (adaptive facets only)
    if (IFACAD != 0) {
        AM_ARTPER(
            MA[LB("KT_IARPER")], MA[LB("KR_RCOURB")],
            MA[LB("KR_DISTAR")], NEDGES
        );
    }
    
    // =====================================================================
    // Assemble curvature radii on facets
    // =====================================================================
    MemoryManager::allocate("KR_RZON12", NFA, 3);
    MemoryManager::allocate("KR_RZON23", NFA, 3);
    MemoryManager::allocate("KR_RZON34", NFA, 3);
    MemoryManager::allocate("KR_RZON41", NFA, 3);
    MemoryManager::allocate("KR_DZON12", NFA, 3);
    MemoryManager::allocate("KR_DZON23", NFA, 3);
    MemoryManager::allocate("KR_DZON34", NFA, 3);
    MemoryManager::allocate("KR_DZON41", NFA, 3);
    
    AM_ASRAY(
        MA[LB("KT_IART12")], MA[LB("KT_IART23")],
        MA[LB("KT_IART34")], MA[LB("KT_IART41")],
        MA[LB("KR_RCOURB")], MA[LB("KR_DISTAR")],
        MA[LB("KR_RZON12")], MA[LB("KR_RZON23")],
        MA[LB("KR_RZON34")], MA[LB("KR_RZON41")],
        MA[LB("KR_DZON12")], MA[LB("KR_DZON23")],
        MA[LB("KR_DZON34")], MA[LB("KR_DZON41")],
        NFA
    );
    
    // Free remaining temporary arrays
    MemoryManager::free("KT_IART12");
    MemoryManager::free("KT_IART23");
    MemoryManager::free("KT_IART34");
    MemoryManager::free("KT_IART41");
    MemoryManager::free("KT_NELART");
    MemoryManager::free("KT_INOAR1");
    MemoryManager::free("KT_INOAR2");
    MemoryManager::free("KT_IARPER");
}